home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / c / jpl_c.zip / POW.C < prev    next >
Text File  |  1986-05-18  |  4KB  |  129 lines
  1. /* 1.0  04-27-84 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1984        *
  6.  ************************************************************************
  7.  *    Programmmed using the algorithm given in:
  8.  *
  9.  *    Coty, William J., Jr., and Waite, William, SOFTWARE MANUAL FOR
  10.  *    THE ELEMENTARY FUNCTIONS, Prentice-Hall Series in Computational
  11.  *    Mathematics, Prentice-Hall, Inc., Inglewood Cliffs, NJ, 1980,
  12.  *    pp. 84-124.
  13.  *
  14.  *----------------------------------------------------------------------*/
  15.  
  16. #include "defs.h"
  17. #include "stdtyp.h"
  18. #include "errno.h"
  19. #include "mathtyp.h"
  20. #include "mathcons.h"
  21.  
  22. /*----------------------------------------------------------------------*/
  23.  
  24. #define reduce(v) ldexp(fint(ldexp(v,4)),-4)
  25. #define K      0.44269504088896340736
  26.  
  27. #define P1    0.8333333333412136e-1
  28. #define P2    0.1249999796500608e-1
  29. #define P3    0.2233824352815418e-2
  30. #define P(a)    (((P3*a+P2)*a+P1)*a)
  31.  
  32. #define Q1    0.693147180559937815e0
  33. #define Q2    0.240226506956777522e0
  34. #define Q3    0.555041084247568661e-1
  35. #define Q4    0.961811769138724104e-2
  36. #define Q5    0.133308101134082075e-2
  37. #define Q6    0.150774061788142382e-3
  38. #define Q(a)    ((((((Q6*a+Q5)*a+Q4)*a+Q3)*a+Q2)*a+Q1)*a)
  39.  
  40. /*\p--------------------------------------------------------------------*/
  41.  
  42. LOCAL double A1[17] =        /* Converted from octal tables    */
  43. {    1.0,            /* in reference above by reading*/
  44.     9.576032806985708e-01,    /* octal values to 45B, and    */
  45.     9.170040432046704e-01,    /* converting to decimal to give*/
  46.     8.781260801866466e-01,    /* A1.  Then this decimal string*/
  47.     8.408964152537131e-01,    /* was converted back to octal.    */
  48.     8.052451659746254e-01,    /* Error in this cycle was    */
  49.     7.711054127039674e-01,    /* added into the converted    */
  50.     7.384130729697489e-01,    /* octal value of the remaining    */
  51.     7.071067811865461e-01,    /* digits, giving A2.        */
  52.     6.771277734684453e-01,
  53.     6.484197773255040e-01,
  54.     6.209289060367418e-01,
  55.     5.946035575013582e-01,
  56.     5.693943173783431e-01,
  57.     5.452538663326258e-01,
  58.     5.221368912137052e-01,
  59.     5.000000000000000e-01
  60. };
  61. LOCAL double A2[8] =
  62. {    2.847357921552306e-15,
  63.     3.151180748043254e-15,
  64.     1.760954860069324e-15,
  65.     7.597361444532664e-16,
  66.     1.065461039677172e-15,
  67.     2.314569550762763e-16,
  68.     2.736854898096527e-15,
  69.     1.721971773273086e-15
  70. };
  71. /*\p*********************************************************************/
  72.     double
  73. pow(x, y)        /* return power function x^y             */
  74.  
  75. /*----------------------------------------------------------------------*/
  76. double x, y;
  77. {
  78.     double c, g, r, xm, u1, u2, v, w, w1, w2, y1, y2, z;
  79.     FAST int iw, p;
  80.     int m;
  81.  
  82.     if (x <= 0.0)
  83.     {    if (x < 0.0 OR x IS 0.0 AND y <= 0.0)
  84.         {    errno = EDOM;
  85.             return -INFINITY;
  86.         }
  87.         else return 0.0;
  88.     }
  89.     g = frexp(x, &m);
  90.     xm = m;
  91.     p = 0;    /* The Ai[] tables are indexed one less than in reference. */
  92.     if (g <= A1[8])
  93.         p = 8;
  94.     if (g <= A1[p + 4])
  95.         p = p + 4;
  96.     if (g <= A1[p + 2])
  97.         p = p + 2;
  98.     c = A1[++p];    /* p is now correct, index one less than ref.    */
  99.     z = ((g - c) - A2[p / 2]) / (ldexp(g, -1) + ldexp(c, -1));
  100.     v = z * z;
  101.     r = P(v) * z;
  102.     r += K * r;
  103.     u2 = (r + z * K) + z;
  104.     u1 = ldexp(ldexp(xm, 4) - p, -4);
  105.     y2 = y - (y1 = reduce(y));
  106.     w1 = reduce(w = u2 * y + u1 * y2);
  107.     w2 = w - w1;
  108.     w1 = reduce(w = w1 + u1 * y1);
  109.     w2 += (w - w1);
  110.     iw = (int) ldexp(w1 + (w = reduce(w2)), 4);
  111.     w2 -= w;
  112.     if (iw > BIGX)
  113.     {    errno = ERANGE;
  114.         return INFINITY;
  115.     }
  116.     if (iw < SMALLX)
  117.     {    errno = ERANGE;
  118.         return 0.0;
  119.     }
  120.     if (w2 > 0)
  121.     {    iw++;
  122.         w2 -= 0.0625;
  123.     }
  124.     m = iw / 16 + (iw < 0 ? 0 : 1);
  125.     p = 16 * m - iw;
  126.     z = Q(w2);
  127.     return ldexp(A1[p] + A1[p] * z, m); /* indices one less than ref. */
  128. }
  129.